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Abstract. In the first course of general relativity, the Schwarzschild spacetime is 
the most discussed analytic solution to Einstein's field equations. Unfortunately, 
there is rarely enough time to study the optical consequences of the bending 
of light for some advanced examples. In this paper, we present how the visual 
appearance of a thin disc around a Schwarzschild black hole can be determined 
interactively by means of an analytic solution to the geodesic equation processed 
on current high performance graphical processing units. This approach can, in 
principle, be customized for any other thin disc in a spacetime with geodesies 
given in closed form. The interactive visualization discussed here can be used 
either in a first course of general relativity for demonstration purposes only or 
as a thesis for an enthusiastic student in an advanced course with some basic 
knowledge of OpenGL and a programming language. 
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1. Introduction 

The visual appearance of a thin accretion disc around a black hole is usually 
determined by direct integration of the geodesic equation from the observer to the 
disc. Unfortunately, this straightforward approach is numerically expensive and an 
interactive exploration of visual effects is cumbersome. However, as light rays are 
independent of each other, the integration can be easily parallelized using many core 
systems, compute clusters, or one of today's graphical processing units (GPUs). If, 
additionally, there is an analytic solution to the geodesic equation, an interactive 
visualization becomes feasible. 

In this paper we take advantage of the analytic solution for light rays in the 
Schwarzschild spacetime together with the highly parallel structure of a GPU for the 
interactive visualization of the optical appearance of a thin accretion disc. 

The first visualization of a thin disc around a Schwarzschild black hole was 
done already by Luminet [T] in 1979. He considered the thin disc being composed 
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of idealized particles moving on timelike circular geodesies around the black hole 
and isotropically emitting light based on the radiation profile by Page and Thorne 
(1974) P]. Here, the mass of the thin disc does not influence the spacetime curvature 
and there is no self-illumination. The bolometric photograph of the disc is dominated 
by the Doppler-shift and by the resulting observed flux that transforms with the 
inverse 4th power of the Doppler-factor. Fukue and Yokoyama (1988) [3J extended the 
visualization of Luminet to multiple wavelengths and compared an X-ray or optical 
photograph with a bolometric photograph. They also discussed light curves when 
an object temporarily eclipses the disc. The visualization by Armitage and Reynolds 
(2003) 4 shows the temporal variability of accretion discs around Schwarzschild black 
holes using magnetohydrodynamic simulations. 

If there is an analytic solution to the geodesic equation of the spacetime under 
consideration, the intersection of a light ray with the disc can be determined 
immediately. This so called emitter-observer problem, which is a boundary value 
problem for the connecting geodesies between the emitting point of the disc and the 
observer, was already discussed in the context of Kerr geodesies by Viergutz (1993) 
or Beckwith and Done (2005) [6]. Details about the analytic solution of geodesies in 
the Schwarzschild spacetime can be found in Chandrasckhar [7] or Cadez and Kostic 
(2005) 8j. A related situation to the thin-disc scenery, namely the analytic observation 
of a single timelike circular geodesic around a Schwarzschild black hole, was discussed 
by Miiller 0. 

The structure of this paper is as follows. In Sec. [2] we present the black hole - 
thin disc scenery and show how the intersection calculation can be reduced to the two 
dimensional situation. In Sec. |3j we briefly outline the analytic solution to the null 
geodesic equation. In Sec. |4j we describe the GPU implementation. We finish with a 
brief discussion. 

Our application is based on the Open Graphics Library (OpenGL) and the 
OpenGL Shading Language (GLSL) [10]. The source code of our cross-platform 
implementation using the QT |llj framework is freely available for Linux and Windows 
systems from http://www.vis.uni-stuttgart.de/relativity. There is also a version that 
can be started directly in a web browser if WebGL [T3] is supported. 

2. Black hole thin disc scenery 

The metric of the Schwarzschild black hole spacetime in spherical coordinates — 
(t, r, d, tp) is defined by the line element 

ds 2 = -(l-^)c 2 dt 2 + ^ , +rW, (1) 
V r / 1 — r s /r 

where dfi 2 = d-d 2 + sin 2 d dip 2 is the spherical surface element, r s — 2GM/c 2 is the 
Schwarzschild radius, G is Newton's gravitational constant, M is the mass of the black 
hole, and c is the speed of light, see e.g. Rindler [15] . 

The geometric outline of the black hole - thin disc scenery is shown in Fig.[l] where 
we use pseudo-Cartesian coordinates as global reference system Q = {e x , e y , e 2 }. The 
transformation between the spherical Schwarzschild and pseudo- Cartesian coordinates 
is as usual: x — rsin$cos(/?, y = rsin^sin^, z — rcosz?. 

In Fig. [I] the Schwarzschild black hole is located at the origin. The disc with 
inner and outer radii i? in and i? out , respectively, is tilted towards the observer with 
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Figure 1. For each pixel P of the virtual image plane of an observer located at 
r = r obs , the intersection Q of the corresponding null geodesic and the thin disc 
with inclination t has to be determined. With respect to the observer's camera 
system, the intersection plane is tilted by the angle ui and the geodesic has initial 
angle £. 



normal lying in the (e x ,e z ) plane, n — (cos t, 0, sin t)g . The observer's position is on 
the e x axis a coordinate distance r oba from the black hole apart. 

2.1. Camera system and geometry reduction 

The reference system C of the observer's camera system, 

C = {e dir = -e (r) ,e right = e (¥;) ,e up = -e (l?) } , (2) 
is given with respect to the local reference frame 



e (*) = — / . 9 t , e (r ) = \ 1 -d r , eM — -d$, et v \ — — ; — -d v (3) 

v ' cWl — r s /r V r r r sin 1/ 

at his position. Here, the observer's four- velocity u = e^) , which means that he is 
static with respect to the asymptotic infinity. The viewing direction e dir points to 
the black hole. The camera's virtual image plane is defined by the vertical field of 
view (fov^) in degrees and the image resolution (res^, res„) in pixels with aspect ratio 
(p = res h /res v ), see Fig. [2] 

For each pixel P = (ph,p v ) of the camera's virtual image plane, a past-directed 
null geodesic with initial four-direction k = (— 1, fc/||fc||)e = ( — 1> k d , k r , k u )c with 
respect to the observer's camera system and 



k = 



1 res h J 2 V resj 2 



(4) 



has to be tested for intersection with the disc. The minus sign in the 0-component of 
k indicates that the null geodesic is past-directed, and the components k d , k r , k u are 
normalized. Because in the spherically symmetric Schwarzschild spacetime the path of 
a null geodesic stays in a fixed two-dimensional plane, we can restrict the intersection 
calculation to the two-dimensional situation shown in Fig. [3] 
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Figure 2. The camera's virtual image plane is defined by the image resolution 
(resh, res„) given in pixels and the vertical field of view (fovt,) given in degrees. 
A pixel P is defined by its pixel coordinates £ [C^res^) and p v S [0,res„). The 
origin for the pixel coordinates is the upper left corner. 




Figure 3. A null geodesic with initial angle £ in the plane V intersects the disc 
plane T> in the point Q. The observer is located at r abB = 15r s and the disc radii 
read i? in = 3r s and i? out = 7r a . The intersection line follows from V n T>. 



Here, the plane V of the null geodesic in pseudo-Cartesian coordinates is tilted 
by the angle to around the e x axis, see also Fig. [2] where tan a; = k u /k r , and = e x , 
= cos lj e y + sin uj e z , e' z = — sin uj e y + cos u e z . The intersection between V and 
the plane of the disc, T>, is given by the straight line 

£(s) = s(-k u tan t,k r ,k u )l . (5) 

The angle <j) between the x-axis and £(s) reads 

i k u t&m 

cos0= . (6) 

v/(fc'') 2 + (l+tan 2 i)(fc") 2 

Hence, an intersection point Q has coordinates (rq = \\£(sq)\\, 4>q = 4>) with respect 
to the plane V . While cf> is fixed by the intersection line, we have to determine the 
coordinate distance tq of Q to the black hole and test whether it lies in the range 
[R in ,R out \. For that, we need the initial direction £ of the null geodesic within the 
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intersection plane, 



cos ^ 



(7) 



i\M\-\\r\\ V 1 + ( fcr ) 2 + ( fcu ) 2 
where r = (1, 0, 0)g with respect to C. The distance rg then follows from the analytic 
geodesic as discussed in Sec. [3] 

As the gravitational frequency-shift between particle and observer depends only 
on the relative distances tq and r obs to the black hole, we can go ahead and determine 
1 + z glav = v/(l - r a /r ob .)/(l - r s /r Q ) (see e.g. Wald O). 



2.2. Disc particles 

Like Luminet pQ , we consider the disc as being composed of idealized particles moving 
on timelike circular geodesies and radiating isotropically. To determine the Doppler- 
shift, we need the angle between the particle's direction of motion and the outgoing 
light ray that reaches the observer. The outgoing light direction will be discussed in 
the next section. 

At the intersection point the non-normalized direction fh = n x £ of a particle on 
a timelike circular geodesic around the black hole reads 



(— k r sin l, — fc"(sin t tan i + cos i), k r cos t) 



9 ' 



(8) 

which we have to normalize and project into the rotated system {e^,, e' y , e' 2 }. 
Depending on the distance tq to the black hole, the particle will have four-velocity 



u = C7 



cy/l - r s /r 



m 



r 



1 - -d r + + m'vd, 



(9) 



with 7 = 1 



1 — P 2 and local three-velocity (see e.g. Muller and Boblest (2011) [T5] ) 
1 



y/2(r/r s -l)' 



(10) 



In Eq. (JoJ) , (m' r , m ,f> , m ,lf ) is the spherical coordinate representation of the normalized 
direction m! in the primed system, where we had to replace m! r by m' r \/l — r 8 /r 
to comprise the Schwarzschild coordinates instead of the spherical Minkowski 
coordinates. Please note that the particle's trajectory appears as an ellipse and thus 
can have a radial component in the primed coordinates even though it is on a circular 
orbit in the global reference system. 



2.3. Radiation of the disc 

The bolometric flux of radiation from a disc around a Schwarzschild black hole is given 

by 



F 



F 



(i?-3/2)i? 5 /2 



R - yfzj2 V3 + v/3/2 



(11) 



with F = 3GMM/(87rrf), accretion rate M, and R = r/r s . The maximum F max /F 
9.167- 10~ 4 is located at i? max ps 4.776. The observed bolometric flux F obs is related to 
the flux F arc emitted by the disc via the transformation F oba = F ero / (1 + z) A . Instead of 
the bolometric flux, we can also consider the observed temperature T ob3 = (F oha / 'a) 1 / 4 
which follows from the Stefan-Boltzmann law with constant a. 
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A detailed discussion of this topic is out of the scope of this article. The interested 
reader might consult the underlying papers by Page and Thorne [2J, Luminet pQ, or 
Fukue and Yokoyama [3]. 

3. Analytic null geodesic 

The intersection point Q between the null geodesic and the disc can be determined 
using the Euler-Lagrange formalism (see e.g. RindlerfTSJ) with the Lagrangian 



C = -[l--)c 2 i 2 + -^—+r 2 $ 2 (12) 
V r J 1 — r s /r 

for geodesies in the I? = ir/2 hypersurface of the Schwarzschild spacetime. Here, the 
dot represents the derivative with respect to the affine parameter A, i.e. f = dr/dX. 
From the Euler-Lagrange equations we obtain r 2 + (1 — r s /r)h 2 /r 2 — k 2 /c 2 with the 
two constants of motion k = (l — r s /r)c 2 t and h = r 2 (j>. For a null geodesic with initial 
angle £ with respect to the observer's reference frame, Eq. ([3]), these constants read 
k = ±cyfl — r s /r obs and h = r oba sin£. The coordinate transformation £ = r s /r of the 
r Euler-Lagrange equation together with C — for null geodesies yields the orbital 
equation 

s = f§V s= ° a -( 1 -0< a with a2 = Cl^|f. (13) 

02 \d(p J sm z £ 
The Mobius transformation £ = (CiT 2 — (2)/(t 2 — 1) brings Eq. ( 13 ) into the standard 



form of an elliptic integral of the first kind, 

^ = +— g— / — - dT (14) 

see e.g. Lawden |16j . The corresponding primitive reads 

(15) 




with module m 2 = (d — C3 ) / ( C2 — C3) an d the three roots Q of the cubic equation 



a — ( 1 — C)C — of Eq. ( 13 ) which can be solved using Cardano's method. Depending 
on the algebraic sign of the discriminant D = a 2 (a 2 /A — 1/27) of the cubic equation, 
these roots read 

Ci,2=p('cos|±\^sm|^ +i Ca = -2 P cos | + i D < 0, (16) 

Ci,2=p('cosh|±iV3sinh|^ +i C3 = -2/)cosh| + i i?>0 (17) 

with the auxiliary variables q = a 2 /2— 1/27, p = sign(q)/3 and cos-0 = q/p 3 for D < 
and cosh-0 = (7/p 3 for D > 0. The order of the roots can be changed at will. 

Since we already know the angle <f> from the intersection line, we need the inverse 



of Eq. (15 1. With the Jacobi-sn function being the inverse of the Elliptic integral 



function T and 4>q = <j)(C — Cobs) from Eq. (151, we obtain 



= C2~Cisn 2 (^yC2^C^>±0o),m) 
l-sn2(iVCW3(0±0 o ),m) ' 
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Now, we can determine rQ = r s /^Q from (18 1 with <f> = <pQ, £ obs = r s /r abs , and £ 
from Q. For that, we have to restate Eq. (18) such that all expressions become real 
using transformation rules as described in Lawden |16j . or we have to do complex 
arithmetics. (Details can be found in the source code.) 

To determine the Doppler-shift, 1 + zd, due to the motion of the disc particle, 
we need the null direction k = k^d^ of the outgoing light ray at (rQ, 4>q). In terms of 
spherical coordinates in the plane of the geodesic, V, the null direction reads 

k = idt + rd r 



1 WS-(i-?)5a r+ ^ as) 



c 2 (l — r s /r) V c 2 V r ' r 2 r 2 

with r — rQ. Then, the Doppler-shift follows from the four- velocity u = u^d^, Eq. |9]), 
the null direction, and the Schwarzschild metric g^ u , 

l + z D = 9^u v . (20) 

Together with the gravitational shift 1 + z grav , we obtain the overall frequency shift 
l + z= (1 + z Krav )(i + zd). 



4. Interactive visualization 



The visualization of the black hole - thin disc scenery, Fig. [Tj is well suited for 
implementation on current programmable graphics hardware. Because light rays can 
be handled independently, we can take advantage of the highly parallel GPU SIMD 
(single instruction multiple data) architecture which results in an enormous speed- 
up compared to traditional serial CPU implementations. At the same time, we can 
combine the simulation and the visualization for an interactive exploration of the 
scenery. By means of GPU programming interfaces like OpenGL, GLSL, CUDA, 
or OpenCL, the visualization algorithm as described in the previous sections can be 
implemented in nearly the same fashion as usual C code with only a minor overhead 
for the interface handling. Here, we use the graphics library OpenGL together 
with the shading language GLSL. The graphical user interface is based on the Qt 
framework [TT], but any other framework with an OpenGL/GLSL binding could also 
be used. 

The starting point is a simple window-filling rectangle (quad) of size res?! x res„ 
that represents the observer's virtual image plane, see Fig. [2j The so called fragment 
program on the GPU then generates a light ray for each pixel (fragment) of the 
virtual image plane, determines the intersection between each light ray and the disc, 
calculates the corresponding frequency shift, and transforms the bolometric flux or 
the temperature. In the last step, either the bolometric flux is mapped to a gray 
value, or the apparent temperature is color-mapped according to the blackbody color 
representation by Wyszecki and Stiles [T71 [TH] . 

Figure [4] shows the visual appearance of the bolometric flux of a narrow thin disc. 
Beside the primary image of the disc, there are also higher order images (lower ring 
and upper thin ring). The lower ring shows the rear bottom side of the disc but, 
because of the rotational symmetry of the Schwarzschild spacetime about the optical 
axis, the image is mirrored. Hence, light rays from the left part originate from the 
right hand side of the disc and are emitted in the direction of motion of the disc's 
particles. That is why the lower ring appears brighter on the left hand side. 
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Figure 4. Visual appearance of the bolometric flux of a thin disc around a 
Schwarzschild black hole. The observer is located at r ob3 /r s = 120 and has 
inclination i = 84.5° to the normal of the disc that has inner and outer radii 
Rm/r s = 3 and iJ ou t/r s = 15, respectively. The camera's vertical field of view is 
fov„ = 6°. For a better printout, we use a gamma correction with 7 = 2. 



5. Discussion 

The main difficulty when using GPU programming is that there are still only 
elementary functions like square root or trigonometric functions available. More 
complex functions must be reimplemented by oneself. In our case, we had to implement 
the arithmetic-geometric mean algorithm by Abramowitz and Stegun [19j for the 



Jacobi elliptic functions of Eq. (18) and some basic arithmetic operations for complex 
numbers. 

Since our code is based on the standard graphics pipeline (OpenGL version 2.0) 
using only the programmable fragment unit (GLSL version 1.10), we are limited to 
single precision floating point calculations. However, for an observer who is not too 
far from the disc, r ob3 < 500r s , this causes no significant artifacts. If double precision 
is essential, the presented algorithm could be easily implemented using CUDA or 
OpenCL. 

We tested our implementation with two different graphics boards. The NVidia 
Geforce 8400M GT reaches about 40 frames per second at a resolution of 600 x 500 
pixels. The more recent NVidia GeForce GTX 480, on the other hand, reaches more 
than 400 frames per second at a resolution of 1000 x 1000 pixels. 

Extensive parameter studies about relativistic scenarios that can be cut down to 
the emitter-observer problem like the black hole - thin disc scenario discussed here, 
definitely benefit from a GPU based implementation because of the highly parallel 
architecture of a GPU. Even more relevant, however, is the usage of the (semi-) 
analytic solution to the geodesic equation that avoids the numerical integration of 
the geodesic equation as well as the time expensive intersection calculation. 
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